clear all; clc;close all;
%Updated: 24/05/2018

a   =load('agrid2.txt');
b   =load('bgrid2.txt');

fm0     =load('fm_ss.txt');
def0=load('gdef1.txt');

na=size(a,1);
nk=1;
ns=na/nk;
nb=size(b,1);

lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

for i=1:na;
    for j=1:nb;
        def(i,j)=def0((i-1)*nb + j);   
    end;
end;


for i=1:nb;
for j=1:na;
	fm(j,i)    =fm0((i-1)*na + j);        
	fm1(j,i)   =def(j,i)*fm(j,i)+(1-def(j,i))*(-2);        
end;
end;


b =-b/(lambda+r);  %change of sign
gdpss=4.8419;
a=exp(a);
b=100*b/gdpss;

figure;
contourf(b,a,fm1,5);
xlabel('debt (% GDP^{ss})');AX1 = gca;box on;
ylabel('y^T');set(AX1,'YTick',-0.8:0.05:1.2);axis(AX1, [-10 40 0.88 1.128]);
